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We calculate the long-wavelength static screening properties of both neutral and doped graphene 
in the framework of density-functional theory. We use a plane-wave approach with periodic im¬ 
ages in the third dimension and truncate the Coulomb interactions to eliminate spurious interlayer 
screening. We carefully address the issue of extracting two dimensional dielectric properties from 
simulated three-dimensional potentials. We compare this method with analytical expressions de¬ 
rived for two dimensional massless Dirac fermions in the random phase approximation. We evaluate 
the contributions of the deviation from conical bands, exchange-correlation and local-fields. For 
momenta smaller than twice the Fermi wavevector, the static screening of graphene within the 
density-functional perturbative approach agrees with the results for conical bands within random 
phase approximation and neglecting local fields. For larger momenta, we find that the analyti¬ 
cal model underestimates the static dielectric function by « 10%, mainly due to the conical band 
approximation. 


I. INTRODUCTION 

The electronic properties of two-dimensional (2D) ma¬ 
terials have been intensively studied in the past decade. 
They offer the opportunity to probe exciting new low¬ 
dimensional physics, as well as promising prospects in 
electronic device applications. Among those interesting 
properties is electronic transport. In the context of elec¬ 
tronic transport in graphene, screening is crucial for elec¬ 
tron scattering by charged impurities^ - — . electron-phonon 
coupling^ - —, or electron-electron interactions^. 

Dimensionality is well known to be essential in deter¬ 
mining the physical properties of materials. Correctly 
describing the physics of 2D materials requires careful 
modeling and definition of the relevant physical quan¬ 
tities. This is particularly true for ab initio calcula¬ 
tions based on plane-wave basis set, as they rely on peri¬ 
odic boundary conditions along the three dimensions. In 
this framework, when simulating low-dimensional mate¬ 
rials, periodic images of the system are necessarily in¬ 
cluded in the calculation. For some physical proper¬ 
ties, the interactions between the periodic images are 
sufficiently suppressed by imposing large distances be¬ 
tween them^i. However, if the electronic density is per¬ 
turbed at small wavevector, long-range Coulomb inter¬ 
actions between electrons from different periodic images 
persist even for very large distances, leading to some 
spurious screening. On the other hand, ab initio cal¬ 
culations have the advantage of describing a complete 
band structure and accounting for local fields. Local 
fields designate electronic density perturbations at wave¬ 
lengths smaller than the unit-cell dimensions^ - — . Ac¬ 
counting for local fields usually requires heavier analyt¬ 
ical and computational worbi^r—. They have been es¬ 
timated in various semi-conductors using first-principles 
calculations^ - — and usually renormalize the screening 


by a few tens of percent. 

The static dielectric function of graphene has been de¬ 
rived analytically within a 2D Dirac cone models — in 
the random phase approximation (RPA). In those deriva¬ 
tions, the role of higher energy electronic states, the de¬ 
viation from conical bands and the so-called local fields 
were neglected. Later, quasiparticle self-consistent GW 
calculations^ of the screening of point charges in neu¬ 
tral graphene seemed to indicate a significant contribu¬ 
tions from the local fields. The general behavior of the 
static dielectric function was found to be quite different 
from the analytical RPA derivation. However, Coulomb 
interactions between periodic images were not disabled. 
There has been some propositions^^, within density- 
functional theory, to correct the contributions from the 
periodic images. More simply, complete suppression of 
those spurious interactions can be achieved by cutting off 
the Coulomb interactions between periodic images^ - —. 
In a recent study of the energy loss function of neu¬ 
tral isolated graphene^, the use of a truncated Coulomb 
interaction (Coulomb cutoff) was implemented in the 
framework of time-dependent density-functional theory. 
It was found that the dynamical screening properties of 
graphene were strongly affected by the spurious interac¬ 
tions between periodic images. 

In this work, we focus on the long-wavelength and 
static screening properties of both neutral and doped 
graphene. We use density-functional perturbation the¬ 
ory (DFPT) as it includes the complete band structure of 
graphene and the effects of local field s 35 ! 36 and exchange 
correlation in the local density approximation (LDA). We 
implement the Coulomb cutoff technique and carefully 
address the issue of extracting two dimensional dielectric 
properties from simulated three-dimensional potentials. 
We then compare our DFPT calculations with the ana¬ 
lytical derivations for the two dimensional massless Dirac 
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fermions within RPA. 

In Sec. IH we set the general background of this work 
by defining the static dielectric function in different di¬ 
mensionality frameworks. In Sec. EU we present dif¬ 
ferent methods to calculate the static dielectric function 
of graphene. This includes analytical derivations previ¬ 
ously developed^— and a self-consistent solution im¬ 
plemented in the the phonon package of the Quantum 
ESPRESSO (QE) distribution. In Sec. [IV] those meth¬ 
ods are applied to both doped and neutral graphene and 
the results are compared. 


II. STATIC DIELECTRIC FUNCTION 

In this section we introduce the quantities of interest 
in the formulation of the static dielectric response. We 
use the density-functional framework within LDA and 
atomic units to be consistent with the following ab initio 
study. Both the unperturbed system and its response to 
a perturbative potential are described within this frame¬ 
work. We start with a quick description of the unper¬ 
turbed system. Since we are interested only in the static 
limit here, we consider a time-independent Kohn-Sham 
(KS) potential^ VKs(r), where r = ( x,y,z ) is a space 
variable. This potential is the sum of three potentials : 

VksM = Kxt(r) + Vh(t) + Vxc(r). (1) 

In the unperturbed system, the external potential V ex t is 
simply the potential generated by the ions of the lattice. 
The remaining potentials are functionals of the electronic 
density. The Hartree potential Vh reads : 

Vv(r) = e 2 J dr' (2) 

and Vxc is the exchange-correlation potential. Since 
the KS potential determines the solution for the density 
which in turn generates part of the KS potential, this ap¬ 
proach leads to a self-consistent problem. When solved 
for the system at equilibrium with no perturbation, the 
ground-state density no(r) is found. 

We now proceed to the description of this system 
within perturbation theory. An external perturbing po¬ 
tential (W ext is applied. This triggers a perturbation 
of the electronic density such that the total density is 
no + Sn , where Sn is the first order response to the per¬ 
turbing potential. Likewise, all the previously introduced 
potentials can be separated in a equilibrium and per¬ 
turbed part. The screened perturbation <5 Vks felt by an 
individual electron is the sum of the bare external per¬ 
turbation SV ext and the screening potential SVu + SV X c 
induced by the density response 5n: 

<W K s (r) = <5I4xt 0) + SV h (r) + SV X c (r). (3) 

This leads to an other self-consistent system^ solved by 
the density response Sn. From this response we can ex¬ 
tract the quantities characterizing the screening proper¬ 
ties of a material. The induced electron density Sn can be 


seen as independent electrons responding to the effective 
perturbative potential <5 Vks : 

Sn(r) = Jdr' X °(r,r')SV KS (r'), (4) 

thus defining the independent particle static susceptibil¬ 
ity x°. It can also be seen as interacting electrons re¬ 
sponding to the bare external perturbative potential: 

Sn( r) = J dr'x(i\r')SV ext {r'), (5) 

thus defining the interacting particle susceptibility 
We can now proceed to further description of the 
screening properties of the material. The static dielec¬ 
tric function is first defined in three- and two-dimensional 
frameworks in order to highlight and clarify their differ¬ 
ences. We then treat the intermediary cases of a 2D- 
periodic system of finite thickness and a periodically re¬ 
peated 2D system, particularly relevant for ab-initio cal¬ 
culations. For those cases, we will determine the condi¬ 
tions in which it is suitable to define a 2D static dielectric 
function. 


A. Three-dimensional materials 

In a periodic system, it is more convenient to work with 
the Fourier transform of Eq. [I] Considering a periodic 
external potential <5I4 xt (r) = <5Kxt(q)e* qr of wavevector 
q, we have in linear response theory: 

Sn( q + G) = E X°(q, G , G 7 )(5Vks(q + G'). (6) 

G' 


Here, reciprocal lattice wavevectors G, G' were intro¬ 
duced. Even though <5U ext (r) only has a q component, the 
electronic density response can include larger wavevec¬ 
tors q + G. Consequently, the induced and total po¬ 
tentials can also have q + G components. Those small 
wavelength components (smaller than the lattice peri¬ 
odicity) in the response of the electrons are called local 
fields^. In a three-dimensional framework, the Fourier 
components of the induced Hartree potential are : 


<W H (q+G) 

^ 3D (q + G) 


A 3D (q + G) 

Ato 

47re 2 


<5n(q + G) 


(7) 

( 8 ) 


where Ug D (q + G) is the q+G component of the Fourier 
transform of the 3D Coulomb interaction. The static 
dielectric constant kq renormalizes the Coulomb interac¬ 
tion depending on the dielectric environment. We focus 
here on an isolated graphene layer, so that ko = 1. This 
constant can also be used in a simple Dirac cone model to 
include the effects of other bands^ic—, though no definite 
value has been proposed. This will be discussed in Sec. 
IIVI Until then, we set kq = 1. The Fourier components 
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of the XC potential are written <5Vxc(q + G). From Eqs. 
0 M and 0 we can write : 


6U K s(q + G) = (5T4xt(q) <5g,o + £U X c(q + G) (9) 
+v 3 c D (q + G) ]T x°(q, G, G')<5U KS (q + G'), 


where <5g.o represents Kronecker’s delta. 

The inverse screening function is defined as the ratio 
of the G = 0 component of the KS potential (the coarse¬ 
grained effective potential) over the external potential: 


e 3 D 


(q) 


^VKs(q) 
<5t4xt(q) ’ 


( 10 ) 


B. 2D materials 

We now wish to work with 2D electronic densities 
Sh(r p ), defined in the {x,y} plane as follows: 

/ +oo 

Sn(r p , z)dz, (11) 

-OO 

where r p is the in-plane component of r, z is the out-of¬ 
plane component. 

We first consider the system usually studied in ana¬ 
lytical derivations, which will be called the strictly 2D 
framework. By strictly 2D, we mean that the electronic 
density can be written as follows: 

6n(r p , z) = Sh{r p )S(z), (12) 

where 5(z) is the Dirac delta distribution. There is no 
periodicity in the out-of plane direction. Considering an 
external potential (il4 x t(q p ) with an in-plane wavevector 
q p , we can define the Fourier transform of the 2D elec¬ 
tronic density Sn( q p + G p ) where G p is a 2D reciprocal 
lattice vector. The Hartree potential 8Vu{r p ,z) gener¬ 
ated by this infinitely thin electronic distribution is three- 
dimensional. We thus separate in-plane and out-of-plane 
space variables to stress the fact that the induced Hartree 
potential does extend in the out-of-plane (z) direction, in 
contrast with the density. Using Eq. 0 and performing 
an inverse Fourier transform in the out-of-plane direction 
only, we find: 

<W H (q P + G p , z) = | M <lr + GpK^+^K 

IQ p i 

(13) 

For our purpose, only the value of the Hartree potential 
where the electrons he <5Vn(qp + G p , z = 0) is of interest. 
Similarly, the KS potential also extends in the out-of- 
plane direction, but we consider only the z = 0 value. 
We can work in a 2D reciprocal space with <5h(q p + G p ) 
and the following potentials: 

<5^H(qp + Gp) = Wn(qp + Gp, z = 0) (14) 

<5Vks(qp + Gp) = <W K s(q P + G p ,z = 0) (15) 

6Uxc(qp + Gp) = <5Vxc(qp + Gp, z = 0) (16) 

<5Kxt(q P ) = ^Uext(qp, 2! = 0). (17) 


Note that since q p is in-plane, <5Kxt(q P ) = <5Kxt(q P , z). 
It is then common practice to use the 2D version of Eq. 
0 with the 2D Coulomb interaction v 2D (q p + G p ) (and 
«o = 1) : 


<%(q P + Gp) 

^c D (q P + Gp) 


— v^ D (q p + Gp)<5h(qp + G p ) 
27re 2 

|q p + Gp| 


(18) 

(19) 


We also define a 2D independent particle susceptibility 
as follows : 


5n(qp + Gp) = ^X°(q, G p , G^Uks^p + G' p ). (20) 

G' p 

Working with the 2D quantities defined above, Eq. 0 
becomes : 

<5Vks(qp + Gp) = <5U e xt(q P ) <5g p ,o + ^Uxc(q P + G p ) 
+^ c 2D (qp + Gp) £ x°(q P , Gp, G;)5U KS (qp + G' p )(21) 

g p 

and the definition of the inverse screening function is 
modified as follows: 


£ 2d( < 1p) 


dVks(qp) 
^Wxt (q P ) 


( 22 ) 


We wish to use this definition in an ab initio framework. 
This raises some issues that we address now. 


C. 2D-periodic materials with finite thickness 


In ab initio calculations, the electronic density extends 
also in the out-of-plane direction. In this section we con¬ 
sider the consequences of a finite out-of-plane thickness 
of the electronic density. We consider now an isolated 
layer with an electron density of thickness 2d. The re¬ 
sults of the purely 2D system should be recovered if the 
wavelength of the perturbation is very large compared 
to d. We illustrate this idea by considering an electronic 
density such that : 


/ +oo 

Sn(r p , z)dz = Sh(r p ) 

-OO 

Sn(r p , z) = 0 if |~| > d . 


(23) 


Using Eq. 0 the 2 = 0 value of the Hartree potential is 
then found to be : 


<W H (q P + Gp) = v 2 c D (q p + G p ) 


(24) 


X 



\^+ G p\\ z \Sn(q p + G p ,z)dz. 


From this equation one can easily deduce that the condi¬ 
tion |q p + Gp|d <C 1 is necessary to obtain results equiv¬ 
alent to the strictly 2D system. Since the largest value 
of |Gp | —1 is only a fraction of the lattice parameter, the 
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above condition can only be fulfilled for G p = 0. The 
q p + G p components of the induced perturbation have 
wavelengths comparable or much smaller than d and the 
thickness of the electronic density cannot be ignored. 
However, as long as |q p |d -C 1, the coarse-grained in¬ 
duced potential can be written: 

<5Vk(q p ) « ^ D (q p )<5h(qp). (25) 

Working with reasonably small perturbation wavevec- 
tors, the 0 = 0 value of the coarse-grained induced po¬ 
tential is equivalent to that of the purely 2D system. It 
is then reasonable to use Eq. [21] at G p = 0 and it makes 
sense to define the dielectric function as in Eq. 1551 


D. 2D materials periodically repeated in the third 
dimension 


In ab initio calculations, in addition to the non-zero 
thickness of the simulated electronic density, an other is¬ 
sue arises. Current DFT packages such as QE rely on 
the use of 3D plane waves, requiring the presence of pe¬ 
riodic images of the 2D system in the out-of-plane direc¬ 
tion, separated by a distance c (interlayer distance). For 
many quantities, imposing a large distance between peri¬ 
odic images is sufficient to obtain relevant results for the 
2D system. However, simulating the electronic screening 
of 2D systems correctly is computationally challenging 
due to the long-range character of the Coulomb interac¬ 
tion. As illustrated in Eq. [T3] the Hartree potential in¬ 
duced by a 2D electronic density perturbed at wavevector 
q p goes to zero in the out-of-plane direction on a length 
scale l/|q p |- For the layers (or periodic images) to be 
effectively isolated, they would have to be separated by 
a distance much greater than l/|q p |. The computational 
cost of calculations increasing linearly with interlayer dis¬ 
tance, fulfilling this condition for the wavevectors consid¬ 
ered in the following is extremely challenging. It is thus 
preferable to use an alternative method. 

In order to isolate the layers from one another, the 
long-range Coulomb interaction is cutoff between layers, 
as previously proposed in such context^ - — . We use the 
following definition of the Coulomb interaction in real 
space: 


v c (r p ,z) 


e 2 0(l z - \z\) 

V\ r p\ 2 +z 2 ' 


(26) 


where 9(z) = 1 if z > 0 and 9(z) = 0 if 0 < 0. The cutoff 
distance l z should be small enough that electrons from 
different layers don’t see each other, but large enough 
that electrons within the same layer do. In other words, 
if d is representative of the thickness of the electronic 
density, we need the following inequalities to be true : 


d < l z < c — d. 


(27) 


to cutoff the Coulomb potential midway between the lay¬ 


ers, L = 


The Coulomb interaction is generally used 


in reciprocal space. Setting l z = | and considering an 
external perturbative potential with in-plane wavevector 
<5f4xt(q P ), the Fourier transform of the above Coulomb 
interaction is written as follow s) 32 ! 33 : 


u c (q P T Gp, G^) — 


47re 2 


|q p + Gp| 2 + G 2 
x 1 - e -l q p+ G pl^ cos(G z l z ) , 


(28) 


where G z is the out-of-plane component of the recipro¬ 
cal lattice vector G. In an ab initio framework, the 3D 
Coulomb interaction u 3D should thus be replaced by the 
cutoff Coulomb interaction v c : 


^%(q P + Gp, G z ) — u c (q p + Gp, G z )5n(<\ p + G p , G z ). 

(29) 

Within the DFT LDA framework, the exchange- 
correlation potential is short-range, such that we can ne¬ 
glect interlayer interactions originating from that term. 
When the Coulomb interaction is cutoff and within the 
region z € [— l z ] +l z \, everything happens as if the system 
was isolated, and it can be treated as the 2D-periodic sys¬ 
tem with finite thickness of the previous paragraph. For 
the layer at 0 = 0, and as long as |q p |d <C 1, we can thus 
work with the 0 = 0 values of the potentials and use the 
definition of Eq. [55] for the dielectric function. 


III. STATIC SCREENING PROPERTIES OF 
GRAPHENE 

In this section we present several methods to calculate 
the inverse static dielectric function of graphene. First, 
the derivation of an analytical expression and a semi- 
numerical solution are presented, following Refs. [2ll427L 
Graphene is treated as a strictly 2D material, its elec¬ 
tronic structure is represented by the Dirac cone model, 
the random phase approximation is used and local fields 
are neglected. Then, we present an ab initio method 
based on the phonon package of QE. This second method 
allows one to relax the approximations involved in the 
analytical derivations. 


A. Analytical and semi-numerical solutions 


When the out-of-plane thickness of the electronic den¬ 
sity can be neglected with respect to the wavelength of 
the external potential, we can work in a strictly 2D frame¬ 
work and Eqs. [51] and [55] can be used. In this section, 
two other approximations are used to simplify Eq. 1211 
Namely, we set <IVxc(q P + G p ) = 0 (RPA) and we ne¬ 
glect the local fields, that is, all G p ^ 0 components. Eq. 
[55] then reads : 


e 2 D 


(q P ) = 


i-^4x°(q P )’ 


IupI 


The interlayer distance can be chosen such that c d 
within reasonable computational cost. Then we choose 


(30) 
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where it is understood that X°(q p ) = x°(q p ,0, 0). In 
a model including only tt — n* bands, the independent 
particle susceptibility is written as follows^—: 


x°( q p) = “2 £. d2k H |(k,s|k + q p ,s')| 2 ^ k ^ k+q: 

W s,s' 


/K 


£ k - £ k+q p 


(31) 


The integral is carried out over electronic wavevectors 
k in one valley around Dirac point K, with a factor 
two for valley degeneracy. The indexes s and s' desig¬ 
nate the 7r or 7r* bands. The occupation of the state 
of momentum k in band s is labelled and is the 
corresponding energy. Within the Dirac cone model, 
a linear dispersion is assumed = s hv F |k|, with 
s = — 1 (s = +1) for the 7r (7r*) band, and vf is the 
Fermi velocity. The wavefunctions overlap is then writ¬ 
ten | (k, s|k + q, s')| 2 = (1T ss' cos($k — $k+q p ))/2, where 
f?k ($k+q p ) is the angle between k (k + q p ) and an ar¬ 
bitrary reference axis. The Dirac cone band structure is 
isotropic and depends only on the norm of the pertur¬ 
bation wavevector |q p |. The numerical implementation 
of this integral in the Dirac cone model will be referred 
to as ” semi-numerical solution”. It has the advantage of 
accounting for temperature effects. In the zero tempera¬ 
ture limit and following the tedious but straightforward 
calculus in Refs. \2 l| [27l . the following analytical forms 
can be found. In the case |q p | < 2 k F : 


where kp = is the Fermi wavevector, if £p is the 
Fermi energy taken from the Dirac point. In the case 
|q p | > 2 k F : 


e2Z>(|q P |) = 1 + 


2e 2 2 k F 
Hvf |q P | 


x 



8 k F 


+ 


IQ P I . — i 

TT^sm 

4/cf 


(33) 



Those expressions are relevant for doped graphene. For 
neutral graphene, we are in the case |q p | > 2 kp, but since 
k F —» 0, Eq. [33] simplifies to : 

e 2C (|q P |) = l + ^. (34) 

The following work aims at investigating the validity of 
those expressions. 


B. DFPT LDA solution 

Several approximations (Dirac cone model, neglecting 
local fields, RPA...) were used in order to derive the pre¬ 
vious analytical expressions. Their validity is not obvious 


in graphene. Ab initio methods like DFPT offer the op¬ 
portunity to relax those approximations^. In this section 
we detail how we obtain the 2D static dielectric function 
as defined in Eq. [55] from DFPT. The issues of the pe¬ 
riodic images and finite thickness in the out-of plane di¬ 
rection are treated as previously discussed. The remain¬ 
ing issues are to apply the adequate perturbation and 
extract relevant 2D quantities. The equilibrium system 
is calculated using the usual DFT plane-waves package. 
At that point, interlayer interactions can be neglected 
in graphene. To study the screening properties, we de¬ 
velop the response of the electronic density to an exter¬ 
nal potential within QE. The code originally calculates 
the induced electronic density in response to a phonon 
perturbation^. Here, we replace the phonon perturba¬ 
tion by the perturbation <5I4xt(q P )- This perturbation is 
constant in the out-of-plane direction and modulated by 
a single wavevector q p in the plane. As shown previously, 
the relevant quantity is the 2 = 0 value of the KS poten¬ 
tial, coarse-grained in the plane <5VKs(q P )- Note that 
the G z ^ 0 components are needed to perform a Fourier 
transform and then take the 2 = 0 value. The number of 
G z elements is limited only by the kinetic energy cutoff. 
We then use the definition of Eq. 1551 


1. Technical details of DFPT calculations 


Our DFT/DFPT calculations were performed using 
the Quantum ESPRESSO distribution^. The electronic 
structure is obtained by DFT calculations within the lo¬ 
cal density approximation^ (LDA). Since the electronic 
structure is calculated without cutoff, it can contain 
some spurious interlayer states above the Dirac point. 
In the calculations, it is thus safer to dope graphene with 
holes to avoid those states. We will assume electron- 
hole symmetry and consider the following results valid for 
both electron and hole doping. We use norm-conserving 
pseudo-potentials with 2s and 2p states in valence and 
cutoff radii of 0.78 A. We use a 0.01 Ry Methfessel- 
Paxton smearing function for the electronic integrations, 
a 65 Ry kinetic energy cutoff, and a 96 x 96 x 1 electron- 
momentum grid. The lattice parameter is a = 2.46 A 
and the distance between graphene and its periodic im¬ 
ages is c = 4.0 x a « 9.8 A. The Coulomb interaction 
is cutoff when calculating the response of the system to 
an external perturbative potential. The results presented 
here were obtained for a perturbation wavevector in the 
direction r —s- K of the Brillouin zone. Identical calcu¬ 
lations were performed in different directions. The vari¬ 
ations on the results were small enough to assume that 
the screening properties of graphene are isotropic. Oc¬ 
casionally, variations from this setup were required and 
will be specified. 
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2. Valididity of the 2D framework 

Now we quickly discuss the validity the 2D treatment 
with respect to the thickness of the electronic density. 
Fig. [T] shows the out-of plane variations of the coarse- 



FIG. 1. (Color online) The induced potential (WH(q P ,.z) in 
the out-of-plane direction at different values of |q p |, expressed 
in units of the distance between the T and K points of the 
Brillouin zone. Calculations were performed for a Fermi level 
of ef = 0.25 eV, taken from the Dirac point. Details of the 
numerical calculations can be found in Sec. mr B II The typ¬ 
ical profile of no(G p = 0,2) is represented. The equilibrium 
density was chosen here to have a common reference for all 
perturbations. 

grained induced potential SVn(qp,z) and the equilib¬ 
rium electronic density no(G p = 0 , z) of a single isolated 
graphene layer in our ab initio framework. We use three 
values of |q p | covering the range of values used in the 
following section. In that range, Fig. Q] shows negligible 
variations of the the induced potential over the extent of 
the electron distribution. The two-dimensional descrip¬ 
tion of the screening properties is thus valid. This range 
of wavevectors covers a large span of situations where 
static screening plays a role. For example, in the case of 
electronic transport we are typically interested in values 
of |q p | on the scale of the Fermi wavevector for relatively 
small doping levels. Thickness effects are negligible in 
this situation. 


IV. RESULTS 


them to the analytical solution (Eqs. I32ll34l labelled ’’An¬ 
alytical”) for the static dielectric function of doped and 
neutral graphene. We identify the contributions of tem¬ 
perature, bands, local fields, and exchange-correlation by 
using different methods. When the analytical derivation 
presented in Sec. IIII Al is used, the Fermi velocity is the 
only parameter needed to define the Dirac cone band 
structure. For consistency with the ab initio methods, 
we use the Fermi velocity obtained in the linear part of 
the DFT band structure, such that Hvf = 5.49 eV-A. It 
is well known that electron-electron interactions increase 
this value by approximately 20% (depending on doping) 
within the GW approximation^. The renormalized value 
is in good agreement with experiments. This renormal¬ 
ization is ignored here, but should be accounted for when 
comparing with experiment. Three intermediary meth¬ 
ods were used to investigate the differences between the 
analytical solution and the self-consistent DFPT LDA 
solution. The first is the semi-numerical method intro¬ 
duced in Sec. IIII Al The independent particle susceptibil¬ 
ity x°(q p ) is obtained by numerical integration of Eq. [31] 
and inserted into Eq. [3D] This solution relies on the same 
approximations as the analytical solution but it can be 
carried out at a chosen temperature (or energy smear¬ 
ing) as long as the integration grid is adequately fine. 
The second is labelled ” RPA” and consists in setting the 
exchange correlation potential to zero within the DFPT 
method. The third is labelled ” RPA no LF” and consists 
in evaluating the DFPT independent particle susceptibil¬ 
ity and inserting it in Eq. [3D] This implies using RPA 
and neglecting local fields, as well as a strictly 2D treat¬ 
ment, since Eq. HU was derived in a strictly 2D frame¬ 
work. This method boils down to the evaluation of Eq. 
E0[ within a more complete ab initio model for the band 
structure. Table [T] summarizes the labels and main char¬ 
acteristics of the various methods used in the following 
plots. 


TABLE I. Summary of the various methods used in the plots 
of Sec. IIVI For each method, we report : i) the treatment of 
electron-electron correlation, LDA refering to the use of the 
XC potential within LDA, ”= 0” meaning that the XC poten¬ 
tial is set to zero ; ii) wether local fields are included (YES) 
or neglected (NO) ; and iii) which band structure model was 
used, the full ab initio band structure or the simpler Dirac 
cone model for 7r — 7r* bands. 


Label 

Exchange-Correlation Local Fields 

Bands 

LDA 

LDA 

YES 

ab initio 

RPA 

= 0 

YES 

ab initio 

RPA no LF 

= 0 

NO 

ab initio 

Analytical 

= 0 

NO 

Dirac cones 


In this section we present the results of the full DFPT 
LDA method (Sec. IIII Bl labelled ”LDA”) and compare 
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FIG. 2. (Color online) DFPT LDA results are plotted with 
and without cutoff of the Coulomb interactions. The in¬ 
verse dielectric function, as defined in Eq. 1221 is plotted as 
a function of the adimensional variable |q p |/|r — K|, where 
|r - K| « i.7 A -1 is the distance between the T and K 
points of the Brillouin zone. The calculations were performed 
for neutral graphene (lower panel) and doped graphene (up¬ 
per panel) with ef = 0.25 eV, measured with respect to the 
Dirac point. In the upper panel, we also represent the scale 
Iqpl/fci? where k,F ~ 0.27|T — K| refers to the Fermi wavevec- 
tor in the doped case. Two interlayer distances were used 
c ss 40 A and c « 9.8 A to convey the dependency of the re¬ 
sults (without cutoff) on that parameter. When the Coulomb 
interaction is cutoff, the results are independent of the inter¬ 
layer distance c. Finally, note that for the neutral case at 
small wavevectors and with cutoff, the results are quite sen¬ 
sitive to energy smearing/grid effects. In this situation, we 
used a 140 x 140 x 1 grid and 0.005 Ry energy smearing to be 
as close as room temperature as manageable. 


A. Importance of cutting off the Coulomb 
interactions 

We begin by presenting the DFPT LDA results and 
pointing out the importance of the Coulomb cutoff in 
Fig. [2] We plot the inverse dielectric function obtained 
with the LDA method with and without cutoff. In the 
latter case, we follow the process of Sec. IIII Bl but the 
original 3D Coulomb interaction vf n is used. Two dif¬ 
ferent interlayer distances are displayed, namely c ss 9.8 
A and c « 40 A . It is clear that interlayer interactions 
play a major role in the screening without cutoff, as a 
strong dependency on the interlayer distance is shown. 
For c w 9.8, the effect of the cutoff is drastic. When the 
interlayer distance is increased, the results without cutoff 
slowly approach the results with cutoff. This is also the 
case in the limit of large wavevector. In general, the re¬ 
sults with and without cutoff are similar when the scale 
on which the induced Hartree potential decreases l/|q p | 
is negligible compared to the interlayer distance c. How¬ 
ever, even using large interlayer distance, the effect of 
cutting off the Coulomb interactions remains significant. 
To obtain accurate ab initio results for an isolated layer, 
it is thus essential to cutoff the Coulomb interactions. 
To give a clearer picture of the effects of the Coulomb 
cutoff, we plot the Hartree potential with and without 
cutoff for two different interlayer distances c in Fig. [3] 
With cutoff, the Hartree potentials corresponding to the 
two interlayer distances coincide exactly with each other 
within the region [— l z ;+l z ], lz being half the smaller in¬ 
terlayer distance here. This confirms that within this 
region, everything happens as if the layers were isolated. 
Without cutoff, in contrast, the Hartree potentials are 
significantly different, stressing the fundamental differ¬ 
ence in the response of systems with different interlayer 
distances. 


B. Comparison of analytical and LDA methods: 
band structure effects 

In Fig. we compare the LDA results (with cutoff) 
to the analytical solution of Eqs. 13211341 The results 
of the two methods are rather close overall. In doped 
graphene, the LDA results are in very good agreement 
(ps 3%) with the analytical method for |q p | < 2 fcp. A 
more pronounced discrepancy (« 10%) is observed for 
|q p | > 2 kp. In the neutral case, a similar « 10% dis¬ 
crepancy occurs for most values of |q p |, but agreement 
seems to be reached in the small |q p | limit. For neutral 
graphene at small wavevectors, smearing plays a signifi¬ 
cant role. Though not plotted here, the semi-numerical 
method is equivalent to the analytical solution when per¬ 
formed with an energy smearing corresponding to room 
temperature. Using the same energy smearing and grid 
as in DFPT to perform the numerical integration of Eq. 
m showed that smearing effects are negligible except 
in the small wavevector limit of the neutral case. For 













FIG. 3. (Color online) The Hartree potentials are plotted 
in the out-of-plane direction with and without cutoff and for 
two different interlayer distances c ~ 40 A and c « 9.8 A 
. The calculations were performed for doped graphene with 
e F = 0.25 eV, at |q p | « 0.32 |T - K| « 1.2 k F - 

DFPT LDA calculations in this regime, we lowered the 
smearing to 0.005 Ry and changed the grid accordingly 
to 140 x 140 x 1 in Figs. [2] and 0] For this smearing, 
agreement between LDA and analytical results is reached 
around |q p | « 0.025|r —K|. Although quite low in terms 
of what is computationally manageable in DFT, this en¬ 
ergy smearing is still large compared to the value corre¬ 
sponding to room temperature. At room temperature, 
we expect that DFPT LDA calculations would show the 
agreement to be reached for smaller |q p |. In the 0 tem¬ 
perature limit, it should be reached for |q p | —► 0. Thus, 
for graphene in general, we can consider that LDA and 
analytical results significantly differ only for |q p | > 2 kF, 
which corresponds to |q p | > 0 in the neutral case. 

To investigate the origin the « 10% discrepancy above 
2 kp, we use the aforementioned ”RPA no LF” method. 
In Fig. |1 this method gives a smaller inverse dielec¬ 
tric constant than both the LDA (« 8%) and analytical 
16%) methods above 2fcp. Comparing the ”RPA no 
LF” and LDA methods indicates that the combined effect 
of RPA, neglecting local fields, and a strictly 2D frame¬ 
work is a « 8% decrease of the results. As mentioned 
before, the band structure model is the only difference 
between the ”RPA no LF” and analytical methods. This 
suggests that the effects of using the Dirac cone approxi¬ 
mation are more sizable (« 16%) but somewhat compen¬ 
sate the other approximations. Overall, we end up with 
the « 10% discrepancy above 2 k F between LDA and an¬ 
alytical method. When setting the exchange correlation 
potential to zero in DFPT, see ”RPA” in Fig|4l the re¬ 
sults are only slightly changed. This means that neglect¬ 
ing the local fields in the plane (what is meant by RPA 
in the derivation of Eq. OiD and out-of-plane (equivalent 


to making the strictly 2D approximation) have more im¬ 
portant effects than exchange-correlation. Although the 
use of an LDA exchange-correlation potential has negligi¬ 
ble consequences for the results presented here, we would 
like to point out that such potentials are derived in the 
framework of a three-dimensional electron gas. Conse¬ 
quently, their relevance in a 2D framework is limited and 
the RPA method might be more reliable than the LDA 
one. 


Iq l/k^ 

F 

0 12 3 4 



Iq l/T-KI 

T 

FIG. 4. Comparison of the static dielectric function of 
graphene obtained within the LDA and analytical methods. 
We use the same axes as in Fig. [2] We also plot the results of 
the ”RPA” and ”RPA no LF” methods. For this last method 
applied to neutral graphene, the point with the smallest |q p | 
was not converged and is not represented. 

A better interpretation of the effects of band struc¬ 
ture can be achieved by comparison of the independent 
particle susceptibility x° from the ” RPA no LF” and an¬ 
alytical methods in Fig. [5] In the |q p | < 2k p regime, 
the screening is dominated by the zeroth-order of x°, 
proportional to the density of states. The linear part of 
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the DFT band structure of graphene is well represented 
by the Dirac cone model. As long as the Fermi level is 
reasonably small (but finite), the densities of states ob¬ 
tained in DFT and analytically are very close. We then 
find a very good agreement with the analytical derivation 
in this regime. In the upper panel of Fig. [5] it is clear 
that a higher-order (in |q p | ) term in x° from DFPT is 
responsible for the gradual disagreement with the ana¬ 
lytical solution as |q p | increases. In the neutral case, the 
zeroth order of x° vanishes with the density of state, and 
X° is always dominated by contributions of higher-order 
terms. For the |q p | > 2ftp regime in general, the first- 
order in |q p | seems to dominate. The susceptibility x° is 
then ruled by interband processes, some of them going 
beyond the range of validity of the Dirac cone model. 

Overall, we find a rather good agreement with the an¬ 
alytical derivation of Refs. I2ll-l2il. This is in strong con¬ 
trast with the conclusions of a previous ab-initio study2& 
of the screening of point charges in neutral graphene. 
Our work differs notably on the use of a Coulomb cutoff, 
and the treatment of ab initio results to extract the 2D 
screening properties of a system that is effectively 3D. 
The authors of Ref. HI state that they checked the neg¬ 
ligibility of the interlayer interactions by looking at the 
effects of interlayer distance on the bands. Such test is 
misleading. Indeed, interlayer interactions are negligible 
on the bands for spacing larger than « 5 A. However, 
as discussed in Sec. EDI the interlayer interactions af¬ 
fect the calculation of the dielectric response when the 
wavelength of the perturbation is comparable with the 
interlayer distance, making the use of a Coulomb cutoff 
essential. We can also comment on the use of the con¬ 
stant /to in Eq. [7] to include the effects of other bands. 
Such a constant is not appropriate since it would affect 
all the orders in x°, including the zeroth order that is cor¬ 
rect. To have an analytic expression quantitatively closer 
to the DFPT LDA results, one should only renormalize 
the contribution from the interband processes. Finally, 
as mentioned before, we used the DFT Fermi velocity in 
this work. One should keep in mind that within the GW 
approximation and consistent with experimental results, 
the Fermi velocity is increased by 20%. This yields very 
similar curves, with a ss 16% increase of the value of 
at large q p , as easily found by plotting the analytical 
expressions. 


V. CONCLUSION 

Definitions of the dielectric function depend on the di¬ 
mensionality. The study of the screening properties of 2D 
materials first requires precise definitions of the relevant 
quantities. After setting such a formalism, we review pre¬ 
vious analytical derivations of the screening properties 
of graphene. We highlight the approximations involved 
in those derivations and propose a DFT-based method 
to overcome them. The DFPT method with Coulomb 
cutoff presented here is general and can be applied to 


Iq l/k 

A p F 
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Iq l/ir-KI 

FIG. 5. Comparison of the independent electron susceptibil¬ 
ity of graphene obtained within the ”RPA no LF” and an¬ 
alytical methods. We use the same axes as in Fig. [2] The 
contribution of intraband and interband processes to x° are 
represented by circles and crosses, respectively. To calculate 
those contributions, we used the semi-numerical solution with 
a small energy smearing (0.001 Ry). The analytical and semi- 
numerical methods are equivalent in that case. 


study the screening properties of other 2D materials. We 
showed that cutting off the Coulomb interactions is es¬ 
sential to recover the screening properties of an isolated 
layer. Our DFPT LDA calculations on graphene lead to 
an inverse dielectric function that is very close to the ana¬ 
lytical form of Refs. f2ll427l for |q p | < 2 kp, and smaller by 
~ 10% for |q p | > 2fcp. Overall, the Dirac cone model in 
a strictly 2D framework, in the zero temperature limit, 
using RPA and neglecting local fields leads to a quite 
accurate and simple analytical expression for the static 
dielectric function of graphene. Smearing effects are neg¬ 
ligible at room temperature and exchange-correlation ef¬ 
fects within LDA are also quite small. Neglecting the 
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local-fields leads to a w 8% underestimation of the in¬ 
verse dielectric function above 2 kp. The largest error 
comes from the Dirac cone model for the band structure. 
This model remains an excellent approximation in the 
|q p | < 2 k,F regime, as long as the Fermi level lies in the 
region where the bands are linear. In the |q p | > 2 kF 
regime, however, the Dirac cone model leads to a k 16% 
overestimation of the inverse dielectric function due to 
the contribution of interband processes probing states be¬ 
yond the Dirac cones. This overestimation compensates 
the local fields effects and the analytical model ends up 
overestimating the DFPT LDA inverse dielectric function 
by « 10% above 2fcp. 
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